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Abstract 



l ^" ^ In this paper, we propose upper and lower error bounding techniques for reduced order mod- 

elhng applied to the computational homogenisation of random composites. The upper bound relies 
on the construction of a reduced model for the stress field. Upon ensuring that the reduced stress 
d satisfies the equilibrium in the finite element sense, the desired bounding property is obtained. 

Ch The lower bound is obtained by defining a hierarchical enriched reduced model for the displace- 

I I ment. We show that the sharpness of both error estimates estimate can be seamlessly controlled 

by adapting the parameters of the corresponding reduced order model. 

^ keywords: Model Order Reduction, Error Estimation, Computational Homogenisation, Proper 

Orthogonal Decomposition, Constitutive Relation Error 
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^ 1 Introduction 

Reduced order modelling is becoming an increasingly popular tool to solve parametrised or time- 
dependent problems {e.g. PHUHl EHl IHl El HH HI] ) ■ Such problems appear in a number of applications 

• • in solid mechanics, including the treatment of uncertainties, structural optimisation and multiscale 

_ ^ modelling. Here, we are particularly interested in the case of nested computational homogenisation 

schemes for random heterogeneous materials {e.g. [40l [TOl [12]) with parametrised micro-structural 
properties. One of the numerical bottlenecks of such approaches is the numerical costs associated 

^ to the statistical volume element problem. Indeed, these costs may be prohibitive if the material 

hetorogeneities are parametrised in order to proceed to their identification or to their optimisation. 

Reduced order modelling proposes to deliver a surrogate model for the solution to a parametrised 
problem, whose evaluation should be inexpensive an accurate. Such techniques consist of two phases. 
A training (or "offline" ) stage and an evaluation (or "online" ) stage. During the training stage, 
the parameter domain is explored, which provides training data that are used to build the surrogate 
model over the parameter domain. In the evaluation stage, the surrogate is evaluated at a particular 
point of interest of the parameter domain. Reduced order modelling techniques differ in the way 
they explore the parameter domain, define and evaluate the surrogate model. One of the classical 
family of reduced order models is based on the response surface method. The solution is evaluated 
at certain points of the parameter domain, and interpolated using explicitly defined shape functions, 
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for instance polynomials. A more advanced class of interpolation techniques such as Kriging, Moving 
Least-Squares or the Reduced Basis Method [3TJ [331 [TT] , do not explicitly define the shape functions 
over the parameter space, but construct them "on-the-fly" during the evaluation stage, based on some 
optimality criteria. Another class of reduced order modelling techniques perform a spectral analysis 
of the training data, and use the part of the spectrum associated to high energy content to build 
the surrogate. This is the case of the methods based on the Singular Value Decomposition (SVD) 
and their extensions (Proper Orthogonal Decomposition (POD) [T3J[5S], Balanced Truncation (see for 
instance Multilinear Singular Value Decomposition, Proper Generalised Decomposition (PGD) 
pn [211 [20]), or on Krylov Subspaces (Moment Matching Method, see for instance [3]). As in the 
case of interpolation-based reduced order models, the surrogate can be defined by an explicit spectral 
expansion over the parameter domain, like in the case of the POD. Otherwise, an implicit definition 
for the surrogate will require some "online" computations to conform to an optimality criterion, like in 
the case of Galerkin-POD [Ml [IHl [T51 [T3] used in this contribution. In any case, the choice of the best 
reduced order modelling method for a particular application strongly depends on the characteristics 
of the underlying parametrised problem. 

Reduced order modelling should be associated to error estimation to control the distance between 
the solution of the parametrised problem (often called "truth" solution) and the solution delivered by 
the surrogate, this requirement becomes even more stringent in the case of nested approximations, for 
instance in the case of multiscale modelling. It is fundamental to understand that there are two types 
of error associated with reduced order modelling: the "offline" and the "online" error. The "offline" 
error is a global distance between the surrogate approximation and the "truth" solution over the whole 
parameter space. Error estimates for this quantity measure the average accuracy of the surrogate over 
the parameter domain. This error can be evaluated by several means, depending on the application, 
which includes cross-validation estimates [T71 [3^1 UHl [S] ) estimates based on spectral analysis of the 
training data (see the review proposed in ^), or more advanced techniques that deliver bounds for a 
particular class of parametrised problems [HI [20] . These estimates are used to control the sampling 
of the parameter domain and the construction of the surrogate. The "online" error is the distance 
between the "truth" solution and the solution delivered by the reduced model at a particular point of 
interest of the parameter domain. The accuracy and reliability of "online" error estimates is crucial 
as they measure the quality of the solution that is actually delivered by the reduced model, which 
is in general not simply related to the average quality of the surrogate]^ In addition, the numerical 
complexity of the "online" estimates should remain low for the evaluation stage to be performed at 
cheap costs. 

This paper focusses on the reliable, accurate and efflcient bounding of the "online" error in the 
context of the Galerkin-POD. In particular, in relation to computational homogenisation, we will 
focus on an elastostatic problem with discontinuous and parametrised elasticity constants, discretised 
by the Finite Element Method. The finite element mesh will be considered sufficiently fine so that 
the discretisation error can be neglected over the reduced order modelling error. The Snapshot-POD 
methodology will then be deployed to extract "offline" an attractive spatial manifold, or reduced space, 
in which any solution to the parametrised problem of elasticity can be accurately represented. In the 
evaluation stage, an optimal solution corresponding to a particular set of elasticity constants can be 
optimally computed by a Galerkin projection of the governing equations in the reduced space. Upper 
bounding techniques for the reduced modelling error have been obtained for such problems in the 
context of the Reduced Basis Method. In [TT], the error estimation relies on the computation of a 
Riesz representation of the parametrised residual, using a fixed bilinear form over the parameter space. 
The bounding property is obtained by weighing the result by a coercivity constant, evaluated "online", 
which is a characteristic of the elliptic operator associated to the parametrised problem of interest. In 
this contribution, we propose to proceed differently by using the concept of the Constitutive Relation 

^One of the fundamental ideas behind the Reduced Basis Method |31| is to hnk the two types of errors by measuring 
the "offline" error using a "max" -type norm in the parameter space. In general, the distinction between "offline" and 
"online" errors is not necessarily relevant if the actual output of the analysis is the quantity used as a target for the 
construction of the surrogate. 
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Error (CRE) , which, in the context of hncar elasticity, only requires to manipulate the concepts 
of displacement and stress admissibilities. 

The CRE is a widely used technique to bound the error associated to a displacement-based approx- 
imation of elasticity problems. In particular, it has been applied to the evaluation of discretisation 
errors in a finite element context |22j , where it coincides in practical implementations with the Equili- 
brated Residual approach (see for instance |H 1351 H] ) . Conceptually, the CRE proposes to construct 
a recovered stress field that is statically admissible, or equilibrated. Applying the constitutive relation 
to the kinematically admissible finite element solution that needs to be verified, one obtains a non- 
equilibrated stress field, called finite element stress field. The distance (in energy norm) between the 
recovered stress field and the finite element stress field is a bound for the error of discretisation. All 
the technical difficulty resides in the construction of the equilibrated stress field [50] . 




Figure 1: Schematic representation of the Galerkin-POD error bounding method based on the Consti- 
tutive Relation Error. 

We extend this idea to the certification of the Galerkin-POD, which will provide a bounding tech- 
nique that is conceptually simple to understand. Here, the reference is the finite element solution. 
Therefore, we first redefine the notion of statical admissibility, and require the recovered stress field to 
verify the equilibrium in the finite element sense. Then, at any point of the parameter domain, we can 
upper bound the error of reduced order modelling by measuring a distance between this recovered field 
and stress field that is directly post-treated from the displacement delivered by the reduced model. 
In order for the recovered stress field to be available at cheap costs in the "online" phase, we build 
a surrogate model for the finite element stress field. This reduced model is completely symmetric to 
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the one developed for the displacement (see figure [T] which can serve as visual guidance for the formal 
developments proposed in this paper). It requires an "online" training from the initial sampling of 
the parameter domain, and an "online" computation to satisfy an optimality condition. While the 
optimality of the reduced displacement is enforced by minimisation of the potential energy, the opti- 
mality of the recovered stress is obtained by minimisation of the complementary energy in the space 
of stresses that are generated by the surrogate model. The technique proposed in this paper is largely 
influenced by, and complementary to, the developments given in '20', where the CRE is applied to 
evaluate the "offline" error of a PGD reduced order modelling technique. 

We will show that the efficiency of the upper bound can be seamlessly controlled by the order of 
the POD expansion of the reduced order model for the stress field. In order to lay the foundation for 
the development of adaptive reduced order models, we also propose a lower bound, which requires the 
"online" solution of a hierarchically enriched reduced model for the displacement and whose effectivity 
can also be controlled. We will also provide some numerical results about the convergence of both 
bounds with the refinement of the initial surrogate. 

The paper is organised as follows. In section[2] we define the parametrised problem of elasticity and 
its discretisation by the Finite Element Method. In section [3) we describe the basics of the Galerkin- 
POD reduced order modelling approach. The error bounding techniques are developed in section 
[4j Finally, the numerical results of the certified reduction technique applied to the homogenisation 
problem will be presented in section [5] 

2 Parametric problem of elasticity solved by the Finite Ele- 
ment method 

2.1 Problem statement: linear elasticity 

We formulate the problem of static equilibrium of a linear elastic structure occupying a bounded 
domain £7 in a physical space of dimension d € {2, 3}. Let M be an arbitrary point of domain fl and 
let X = xi Ci + ... + XdCd be its cartesian decomposition in reference frame TZ = (O^^, Cj^, §2, 63). We 
look for a displacement field u £ U{n) = H^{rt) that satisfies the Dirichlet boundary conditions u — w 
on the part diY" of the domain boundary dil. Any displacement field that satisfies the conditions of 
regularity and the Dirichlet boundary conditions is said to be kinematically admissible and belongs to 
space I4^'^{n) C I4{n). We introduce the Cauchy stress tensor field a which belongs to a space S{iV) of 
suflSciently regular second-order tensor fields. A density of tractions t is applied to the structure on the 
part 9ri* :=i9r2\911^ of the domain. A density of forces denoted by b is applied over il. The principle 
of virtual work expresses the equilibrium of an arbitrary stress field belonging to 5^'^ (£7) C S{il) as 
follows: 

\fu* eu^'^'°{n), - [ g:^{u*)dn+ [ b-u*dn+ [ t-u*dr = o, (i) 

where U^'^'^{il) = {v £ U{n)\v\QQ^ = 0}. In the previous equation, e{u) :— ^{"^u + Vu^) is the 
symmetric part of the displacement gradient. The solution to the problem of elasticity is an admissible 
pair (w, o;) € U^^{^) x iS^'^(il) that verifies the isotropic linear constitutive law 

^ = A Tr(i(u)) + 2 ^i(u) := D : g{u) , (2) 

where A and /i are the Lame elasticity constants, and D is the fourth-order Hooke's elasticity tensor. 
The inverse of this constitutive law reads 

|(M)-^Tr(^)|,-;|^:=C':^, (3) 

with E and 1/ the Young's and Poisson's moduli respectively, which are linked to the Lame constants 
by the relationships A = (i+^^jla,,) and /i = ^l^h)- 
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By substitution of the constitutive law into the principle of virtual work, the problem of elasticity 
can be recast in the primal variational form (displacement approach) 

Find uGU^'^in) such that Vu* eU^'^^°{n), a{u,u*) = l{u*) , (4) 

where the symmetric bilinear form and the linear form associated with the problem of elasticity are 
respectively defined, for any fields v and u* in U{fl), by 

a{v,u*)= [ §^{u*):D:g{v)dn, l{u*) = [ b- u* dV + [ t-u*dr. (5) 
Jq Jn Jdn* 



2.2 Parametrised problem of elasticity 

We consider that the input characterising the problem of elasticity arc functions of a finite set of 
scalar parameters (Mi)i6li.nn] that are ordered in a parameter vector /i e . Let V C M"^' be 
the domain of admissibility of n. More Precisely, the force densities b and t, the Dirichlct boundary 
conditions w and the Young's and Poisson's ratios E and v arc functions of parameter /i. The domain 
fl and its boundary split dfl = 951^ U dfl^ could also be parametrised, but this would lead to additional 
difficulties in terms of reduced order modelling, which will not be addressed in this contribution. 
The field variables are now formally redefined as functions of the parameter: 

u: V ^ U{n) g: p ^ s(n) 

H I— >■ u(At) jJL i-> a:(/x) , ' 

This definition is extended to all parametrised variables and spaces. If not stated explicitly, a symbolic 
expression of the type /(/i) G -^(m) will imply the definition of a function / with inputs in V and 
values in space (J^^p J^(/u). A symbol of the type f{y,z;n) will additionally imply that the values of 
/ arc functions of variables y and z. 

For a given fi, a kinematically admissible displacement field will be sought in W^'^(r2;/x) = {v G 
U{il) \ v\dn-" = w{lA}- Similarly, a statically admissible stress field a{n) G »S^^(ri;/z) will satisfy the 
parametrised principle of virtual work: 



t(/z) ■u*dT = 0, (7) 
an* 



\/u* gU^'^'°{Q), - [ g:g{u*)dn+ [ b{ii)-u*dn+ [ 
Jn Jn Jd 

By constraining the stress fields to satisfy a priori the constitutive equation 

^(m) = -D(m) : i (uilj)) , (8) 
the parametric problem of elasticity can be written in the primal variational form, for any fj, € V: 

Find u{p) e U{n) such that Vm* e U^'^'°{n), , , 

a(M(/f),M*;/f) = l{u*;i^) , 

where the parametrised bilinear and linear forms associated with the problem of elasticity read, for 

any {v,u*) in {U{^))^, 

aiv, u*;p)= / i(M*) : -D(^) : g{v) dQ , l{u*\^= j b{fi) ■ u* d^ + / i(/i) ■ u* dT . (10) 
Jn Jn Jdn^ 

2.3 Quantity of interest and input-output map 

Let us consider a set of scalar quantities of interest ordered in a vector Q{fi) = Q{u{i^)) € M. Our goal 
is to construct a surrogate for map Q. 
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The basic strategy appUed in such cases is to compute the quantity of interest corresponding to 
certain (well-chosen) parameter values in fi G V = {/J,f ... ,Mn,}j where 7^ is a subset of V (called 
training parameter domain or snapshot parameter domain), and interpolate the quantity of interest 
in some ways over the entire parameter domain V. Of course, the method of interpolation should be 
driven by considerations of optimality, stability and quality control. 

In particular, interpolating the solution u and not directly the output over the parameter domain 
is usually preferable. Notably, and this is the point that we focus on in this contribution, reliable error 
estimates can be obtained for a certain class of parametrised boundary value problems. We will use 
the Galerkin-POD approach, which will be described in section |3] 

2.4 Finite element discretisation 

We approximate the solutions to the parametrised problem of elasticity by making use of a classical 
finite element discretisation U^{il,) C U{n) of space U{il) [23 H]. As we assumed that domain Q does 
not depend on parameter /i, we can construct a unique finite element space for all the realisations of 
the parametrised boundary value problem. More precisely, the finite element space will be such that 

U'\n) = {veUifl) I Vj G {1, espan((iV,),eii,„„i)} , (11) 

where Vj denotes the j*'' component of vector field v and functions (-/Vi)ig|i are compactly supported 
finite element shape functions belonging to U{fl). 

Let U^ °{n) :=Z^''(f7) n W^^'°(r2) be the space of finite element fields that vanish on dn^" and let 
uP(/i) be a particular field of U^'^{n;ii), for any fi V. The finite element approximation u^(/i) of 
u(fi) is the solution to the following variational problem: 

Find u^ifj.) e U'^^^in) + {mP(m)} such that Vu* e W'''°(f7), 

a(M^M),M*;M)-;(M*;M). ^ ' 

We will assume in this paper that the parametrised Dirichlet boundary conditions conform to 
the finite element space, which means that U^-^'^{n;fi):^U^{n) nU^'^{n;fi) 7^ {}. In this context, 
uP(^) is naturally chosen in the finite element space U^{il,) and we will use the alternative notation 
Mh'P(^) = vPifj,), where u^'^ifi) S U^^^'^ifl; fj,). 

For any n GP, Problem ( |12[ ) can finally be recast in the form: 

Find e U'^'^i^) such that Vu* E Z^^'"(f7), 

a(M^'°(M),M*;A^) -/(M*;A/)-a(M'''P(M),M*;A^), ^ ' 

and the finite element solution is obtained making use of the lifting identity u^(/i) = m'''°(m) + w'^'P (/i) . 
Notice that the finite element solution obtained in this fashion is kinematically admissible, which is 
important for the remainder of the developments. 

In the following, we assume that the finite element space is sufficiently fine so any measure the 
finite element error e^{n) :=u.(/i) — y^{fJ.) is small for all f.i G V. 

3 Galerkin-Proper Orthogonal Decomposition 

The idea of projection-based reduced order modelling relies on the fact that the solutions to parametrised 
boundary value problems are often found to lie in low-dimensional spatial subspaces engendered by 
the parametric dependence. In the training stage, several methods can be used to capture the at- 
tractive subspace numerically, amongst them the immensely popular Snapshot (or Empirical) Proper 
Orthogonal Decomposition I34j. In the evaluation stage, the governing equations corresponding to 
any parameter of interest are projected in the reduced space obtained "offline" . In the case of elliptic 
parametrised boundary value problems, such a projection by means of the Galcrkin method yields an 
optimally interpolated solution, from which the quantities of interest can be extracted. 
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3.1 Projection-based reduced order modelling 

In order to introduce the approximation of the displacement y}^ using a reduced space, let us first recall 
the lifting of the finite element solution over the parameter space 

y^ieV, y}'i^I)^y}''"i^I) + y}''Pi^i). (14) 



In equation (14), u^'" is a function of /x with values mU^'°{i}), which means that the displacement fields 
vanishes on d^l'^ for any fi G P, while the values of function m'^'P belong to the finite element 
space U^{n) and satisfy the Dirichlet boundary conditions of the parametric boundary value problem 
exactly. We assume for now that the lifting w'^'P is known, and concentrate on the approximation of 
the values of the homogeneous remainder u^''^ using a reduced space. 

Let us introduce a basis (0.)ig|i^„_^l € ^ of the representative subspaceZ^'''°(ri) cU^'^{il.) 

in which any value of the displacement u^'*^ will be approximated. This basis is also assumed to be 
known at this stage. Using these notations, and for any /i G T', we look for an approximation yJ;{^) of 
y^{lJ,) to the parametrised problem of elasticity in the form 

y}'i^I)^^/{^i):=u'•°{^J)+y}'•P{^I), where m'^^m) = X] a^(Ai) • (15) 

i=i 

The functionals (Q!i)ie|i,n^| are interpolation weights, called "reduced variables". The interpolation 
weights can be optimally computed by using a Galerkin formulation of the elasticity problem in the re- 
duced space, for any admissible parameter value. The projected problem corresponding to an arbitrary 
£V reads 

Find u''°iiJ.) e W'°(rj) such that \/u* e W'^ifl), 

a(M''°(Ai)i w*; m) — Km*'i a*) — a(u^'^(/^)i w*; m) ■ 



Substituting (151 into ( |16[ ), and using the explicit form of the bilinear and linear forms of the 
parametrised elasticity problem, one obtains the following (small) system of coupled equations in 
the reduced variables (Q!i(M))jeIi,n^]: 

K'-(Ai)t^(M) = r(/i)+r'P(M), (17) 

where the a(/i) G M"* is the unknown vector of weighting coefficients. The algebraic operators used 
in the previous expression are defined by: 

V(*,j)e [l,n^F, K^,-(M) = a(^^.,0^;M) 

Vjell,n0l, F5(/i) = ;(^^.;m) (18) 

Vj e ll,n4, F;'P(/i) - -a(M'''P(/i),^^.;M) • 



The linear system of equations ( 18 1 can be solved inexpensively in the "online" phase, for any 
parameter value of interest. The result is the approximate field m''(a*) = J2i=i + y^^'^if^), 

from which the approximation of the output (5(u'^(/i)) can be estimated. However, in order to be able 



to compute each of the terms of system ( 18 ) efficiently, some additional assumptions on the form of 
the parametrised boundary value problem are required, which is detailed in the next section. The 
computation of the reduced basis itself, and the definition of the lifting w'^'P is done in the "offiine" 
phase, which will be detailed later on. 

3.2 Case of parametric elasticity problems admitting an afRne form 

It is important to remark that a reduced order modelling technique cannot reach the expected numerical 
efficiency if the complexity of some of the operations performed in the "online" phase depend on the 
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dimension of the underlying finite element space. The usual framework is to deal with problems that 
naturally admit an afhne or radial form (see [11] [50] for instance) . When the problem does not admit 
an affine form (e.g.: moving domains, nonlinear problems), an approximation of the initial governing 
equations must be performed in order to retrieve a reducible problem in affine form (see for instance 
0,32, 28, HEHH]). 

We will suppose in this contribution that the parametrised problem of interest (131 admits a natural 
affine form, which reads mathematically: 

aiv,u*;fj.) = ^aj(w,w*)7f(^) 

"1 

Vmg P, V(m*,i;) g = J2^^iu*)lli^i) (19) 

1=1 

rip 

a(H'>'P(/i),M*;M) :=r'P(li*;M) ^Y.'^'" {2^)11^" it) , 



where (ai)jg|]^ „ | are bilinear forms (not necessarily symmetric nor positive definite) and (7f)jg|2 „ | 
are explicitly known functionals of the coefficients, (L) . „, , are linear forms and ('y]) . „, , are 
explicitly known functionals of the coefficients and „ | linear forms and „ | are 

explicitly known functionals of the coefficients. Such a representation of the parametrised boundary 
value problem is obviously at hand if the parametrised data can be written in the form 







D{x,fi) 


rid 

1=1 






t{x, fl) = 


i=l 


V^i e V, Vx 


e 


b{x, fi) = 


rib 






w{x, fj.) 





(20) 



1=1 



where the notations used are consistent with those introduced for expressions (19). 

In this context, assembling the projected system ( [T8| in the "online" phase can be done efficiently. 
The terms corresponding to each of the summands of the affine expansions ( 19 1 can be pre-computed. 
For instance, assembling the projected stiffness K'^(/i), which is an operation of numerical complexity 
a priori related to the dimension of the finite element space, can be written as follows: 

fc=i (21) 
where, V fc G [1, M, V (*, j) G p, n^f , K', = afe(0,, 0,) / eU):Dk:e{l^)dn. 

Operators (K'^ ) can be precomputed "offiine". In the "online" phase, the assembly is simply 

V— fc/fc6ll,r!dl 

done by computing a linear combination of these operators with coefficients {if if^)) ^^^^i „ |- The same 
technique can be used to assemble F''(/i) and F''P(/i), for any /i of interest. As a consequence, the 
numerical complexity of the "online" phase only depends on nj, nt, rib, and of course on the 
dimension of the reduced space. 



8 



3.3 Non-homogeneous Dirichlet boundary conditions in projection-based 
reduced order modelling 

Applying non-homogeneous Dirichlet boundary conditions in projection-based reduced order modelling 
is not a trivial task. In [] , the authors propose to use the usual boundary lifting performed in a finite 
element context. We propose an alternative approach based a global lifting that is (i) consistent with 
the choice of global basis vectors to define the reduced space and (ii) consistent with the 

construction of the dual reduced space used for error estimation, as will be explained later on. 

Our lifting technique makes use of the previously assumed form of the prescribed displacements 
( pO| ). In the "offline" stage, we compute a set of finite element displacement fields ('0j)ie|i,n„] S 
{U'^in))"'" that satisfy 

This set of fields is obtained by solving riw standard finite element problems "offline" with non- 
homogeneous boundary conditions. The choice of G is arbitrary. 

The finite element lifting function u'^'P, which satisfies the Dirichlet boundary conditions for any 
fj. € V, can now be defined by the affine expansion 



Wfier, = (23) 



3.4 Reduced spaces obtained by the Snapshot Proper Orthogonal Decom- 
position 

The purpose of the Snapshot Proper Orthogonal Decomposition |34j is to deliver an orthonormal 
reduced basis (0j)ig|i,ri^] of given cardinality such that the £^(f2)-projection of a set of values of 
the "truth" solution to the parametrised boundary value problem in the space generated by this basis 
is minimised. In our context, this optimisation problem reads: 



where J 



argmin 



J 



ell,- 



-E 



£2(o) 



611/ 



(24) 



In the previous equation, < . , . >c^{n} denotes the usual scalar product of C^{fl), || . ||£2 denotes the 
associated norm and Sij is the Kronecker delta symbol. The snapshot parameter domain V = {^^ | i G 
|l,ns]} C 7^ is a discrete set of Us parameter values that are chosen in V. The set of corresponding 
values of is the so-called snapshot S° = {v}^''^ (fi) \ /i G V}, which is computed "offline" using the 
"truth" finite element approximation. 

The solution to this classical problem of optimisation can be obtained by solving the eigenvalue 
problem 



|i^7 = A7 where ill, ^ {u'^'" ifA) ^ if^)) c^n) V(*,j)ell,n, 



(25) 



After ordering the eigenvalues (Ai)ig|i_„^| of in descending order and denoting by (7^)ig|i_„j the 
corresponding eigenvectors, the solution to the Snapshot POD optimisation problem (24) is given by 



j=l 



a/a" 



(26) 
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where jij denotes the j component of eigenvector 7.. For a given snapshot 'EP, the cardinal- 
ity ricj, of the reduced basis is usuahy chosen "offline" such that the value of the objective function 

J is small enough, which means that the correlated information contained in the snap- 

shot has been correctly captured. Choosing the snapshot itself is a more difHcult issue. Qualitatively, 
the sampling of the parameter domain should be fine enough to capture the variations of y}^''^ induced 
by the parametric dependence. In practice, this requirement can be validated by empirical means such 
as cross-validation [T71 [U [321 US] i or using more advanced numerical methods [121 133 HO] to estimate 
some measure of the "offline" error over the entire parameter domain. 



4 A posteriori error bounding for projection-based reduced 
order modelling 

For a given /i of interest, the distance between the finite element solution y}^{n) and the displacement 
delivered by the Galerkin-POD needs to be evaluated. We will employ the Constitutive Relation 
Error [23] to bound it from above. This concept requires the availability of a stress field that satisfies 
the discrete principle of virtual work, or equilibrium in the finite element sense. Such a field can be 
obtained by "offline" construction and "online" evaluation of a POD-based surrogate for the finite 
element "truth" stress over the parameter space. This surrogate is similar to the one developed for the 
approximation of the finite element "truth" displacement, which makes the upper bounding technique 
conceptually simple. 



4.1 Definition of the interpolation upper error bound 

We consider a realisation of the parametrised problem of elasticity corresponding to arbitrary an 
parameter fi oiV, which is solved approximately using the reduced order modelling technique described 
in section [sj The reduced order model delivers a kinematically admissible displacement field u^{i^) € 
U^'^'^{p,). However, the stress field (T''(/i) ■.— D{^) : e(u'')(/i) does not satisfy the equilibrium in the 

finite element sense a priori. If it did, the solution to the "truth" parametrised finite element problem 
M^(^) would be at hand. The idea behind the proposed error bounding technique is to post-process 
a so-called "recovered" stress field ct(^) € that is equilibrated in the finite element sense. This 

admissibility condition reads: 

^fu* eU^'^iri), - [ giii) ■■ giu*)dn+ [ b{^i)-u*dn+ ( t{fi) ■ u* dT = . (27) 
Jn Jn Jan* 

We denote by S^-^'^{n;n) the space of stresses satisfying the parametrised equilbrium in the finite 



element sense (27) 



Theorem: The distance between the stress field fT''(/i) G 5(ri) obtained by direct evaluation 

of the reduced order model and the recovered stress field ^{fi) S S^'^^{il,; /i) can be used to bound the 
error of reduced order modelling e''(/i) :=m^(^) — u'^in) as follows: 

- mWcM > l|e'-(M)llD(M) , (28) 



where, ||m*||_d(;^) = [^/o^(^*) • • £{u*)dnj is the energy norm associated to an arbitrary dis- 

1 

placement field u* £ U^''^{fl), and ||a;*||c(/j) = ^ Jq 1* • Qif^) ■ 1* is the energy norm associated 

to an arbitrary stress field a* G 5(0). 
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Proof: The proof is a straightforward extension of the one used to bound the discretisation error in 
the context of finite element approximations [23i . We start by expanding the distance between the 
reduced stress field (f^ifJ,) and the recovered stress field ^(/i) by using the trivial identity 

Wg^if,) - = II (^'-(m) - + m) llc(M) . (29) 



lk'(/i) -|(/i)ii|(^) = M{f^)-u\mm + \\g\f^)-m\\ci,) 

- Km)) : k (^' (a^)) - £ (^^'(a^))) ■ 



where q^{ii) := D{^) : e(u'^(/i)) is the "truth" finite element stress field {i.e.: the stress field which 

would be obtained without reduced order modelling). Then, by using the constitutive relation and the 
definition of the energy norms given previsouly, we can expand identity (|29[) as follows: 



(30) 



Recall that both the finite element stress field and the recovered stress field are equilibrated in the 
finite element sense (equation (27)). Therefore, the following identity holds: 

Vm* €Z^'^'°(r!), / {g^ilj)-gil£)) ■■g{u')dn = 0. (31) 
Jn 

As both the finite element displacement field and the displacement field obtained by solving the reduced 
order model belong to the finite element space L(^{ft) and satisfy the Dirichlet boundary conditions, 
we obtain that the last term in (30) vanishes, which yields 

. 2 



(32) 



The last term of the right-hand side of the previous equation is positive, which concludes the proof. □ 

A first important remark is that the Galerkin orthogonality property a{(f{ii), u*', fi) = for all u* € 
U^''^{n) is not used, which means that the bounding property would hold true if an interpolation method 
other than Galerkin reduced order modelling had been used to interpolate over the parameter domain 
{e.g. Kriging, interpolation over a priori defined polynomial bases, etc.). However, u^{iJ-) would still 
be required to be kinematically admissible. 

A second important remark is that the efficiency of the error estimate depends on the distance 
||(T^(/i) — CT(/i)||(7(^),(i.e.: the distance between the recovered stress field and the "truth" finite element 

stress field), for any fi dV, which needs to be kept in mind when constructing the equilibrated stress 
field. If the recovered stress field is finite element stress field, the CRE estimate j/"p(//) coincides 
with the exact error. Of course, the finite element stress field is not affordable in the "online" phase. 
Therefore, the next question is ''how can wc compute a recovered stress field that is 

• a good approximation of the finite element stress field, 

• equilibrated in the finite element sense, 

• of "online" numerical complexity approximately equal to the numerical complexity of evaluating 
the reduced order model?" 



4.2 Principle of the construction for the equihbrated stress fields 

In order to address the previous question, we propose to build a POD-based reduced order model for 
the finite element stress. The technique is similar to the one used to build the reduced order model for 
the displacement field. First, we use the affine representation of the external load to build a statically 
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admissible stress field over the whole parameter domain. Then, the complementary part of the finite 
element stress field defined over the parameter domain is sampled, and a basis for the subspace spanned 
by these samples is extracted using a singular value decomposition. These operations are performed 
"offline" . In the "online" stage, for a given parameter, the coefficients associated to the reduced basis 
functions are optimally obtained by solving a projected problem. 

Specifically, we start by splitting the finite element stress into two parts, as follows 

e V, g}\ii) = +^'^'P(/^) . (33) 

The first part a^-'^ belongs to a space 5^'°(il) of stress fields satisfying the homogeneous equilibrium 
conditions associated with the "truth" finite element problem, which reads 

V^t e 7^, Vu* e Z^^'"(f7), [ g^-^iii) g{u*) dn = . (34) 

Jn 



The second part of split (331 is a particular stress field CT^'P(/i) G iS^'^'^(f2) that satisfies the equilibrium 



in the finite element sense (equation (27)). (t^'P will be explicitly defined as a function of the parameter. 



while the complementary part a^'^ will be approximated using the Snapshot-POD: 

V/ieP, ^'^(a^)«|(m):=^'''°(a^)+^'''P(m), (35) 

where the approximate stress a^''^{fi) E S^'''\n) is such that it satisfies the homogeneous equilibrium 
equations in the finite element sense for any fi inV. 

It is clear that within this framework, the recovered stress field ^(/i) satisfies the equilibrium in the 



finite element sense (27) over the entire parameter domain. 



4.3 Riesz representation of the parametrised static load 

Let us first construct the particular stress a^'P associated to non-homogeneous equilibration condi- 
tions. We make use of the assumed affine form of the Neumann boundary conditions and body forces 
given in expression ( pO| ). In the "offline" phase, we compute a set of global finite element vectors 
(V'i)ieli, nt-i-nbl G (W^''^ (f2)) Corresponding to the summand of the affine form identities by solv- 
ing successively the finite element problems: 

V^e II, ntl, Vu* eZ^h'°(r!), a(^„w*;^o) = / ■ u* dV 

(36) 

Vze II, nbl,Vu*eZ^'^'*'(r!), a(^,+„,,w*;/io) = / lru*dn 



n 

Notice that we have chosen to enforce homogeneous Dirichlet conditions for this series of finite element 
problems, which ensures a certain regularity of the procedure (the parametrised bilinear form a( . , . ; /io) 

is positive definite over iU^''^{Q?))^)- Parameter /j,o is the one that we also used to define the lifting 
u^'P. This choice has been made for the sake of simplicity. ^ 

In the "online" phase of the error estimation procedure, the fields ['4ii)i^ii^nt+n^l being at hand, 
and for a given parameter fj, G V, we evaluate the field v^'^ifJ.) € U^'^{^1) defined by the formula 

t^^in) = E + E ^^+"* ^^^^H) ■ (37) 

We can then verify that the stress field 

a'^-P(M):=^(/io):£(M''P(A*)), (38) 
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is statically admissible in the finite element sense, for any fi ^V. Therefore, the quantity ||(T'^'P(/i) — 
a' (/i)llc(^) is an upper bound for the error measure ||e''(/i)||£i(^). However, in the case where the body 

load and applied tractions are zero, this bound is trivial of no practical interest. In the general case 
anyway, this bound should be sharpened by computing the complement a;''°(/i) to the recovered stress, 
which is done next by making use use of an additional spectral analysis of the training data. 



4.4 Snapshot POD for the finite element stress 



The "offline" computation of the snapshot delivers a set of admissible stress field in the finite element 
sense S := {ct'^(^) | /i G V}. After subtracting the corresponding values of the previously defined 

statically admissible stress component ct'^'P, we obtain the set 5" := {(T^'''(/i) | /i S V} of stress fields 
satisfying the homogeneous equilibrium equations in the finite element sense. In the "online" phase, 
our purpose is to compute the correction corresponding to an arbitrary parameter ^ E V 

as an optimal combination of these fields in the sense of the minimisation of — CT'^(/i)||c(Ai) = 

||(t'''*'(/x) — o^'^{^J)\\c(p.)■ In other words, we aim at maximising the efficiency of the error estimate. 

In order to do so, we compute a singular value decomposition of the finite element stress fields 
contained in the snapshot 5°. We look for a set of basis tensor fields (</'.)iG|i.fi^] G (5(51))"* that 
are orthogonal with respect to the inner product < . , . >c(aio)~ /n ■ • Q(t^o) ■ ■ dfl of space S{fl), and 
that are solution to the following optimisation problem: 



.)iell,n, 



where J 



argmm 

,6(5(0))"*, <0*,0*>c(^o)=5i, 



J 



((f). 



((f): 



611/ 



-i I C'(aio) =1 



6llv 



C'(mo) 



(39) 



This problem is similar to problem problem (24) exept that we deal with second order tensors, 
and that the optimality of the decomposition is defined in the sense of the weighted norm || . ||c(mo) 
associated to < . , . >c(/jo) h^stead of the usual Frobenius norm. The energy norm used to define the 

optimality of the POD is evaluated in /iq to avoid adding unnecessary parameters to the bounding 
technique. Similarly, we have assumed that the sampling of the parameter domain performed to train 
the surrogates for the stress and and for the displacement are the same in order to keep the methodology 
simple. 



The solution to optimisation problem ( 39 ) is obtained by solving the eigenvalue problem 



H7 = A7 where H^, = (^'^■°(MD:^'''"(M^-))e(^„) V (*, j) G [1, n^l 



(40) 



After arranging the eigenvalues (Ai)ig|i of H in descending order and denoting by (ci ^i<^\\,n^\ the 
corresponding eigenvectors, the solution to the Snapshot POD optimisation problem (39) is given by 



Vzell,n^l, f=E^''°(^' 



7». 



(41) 



J = l \/A,; 

Finally, the recovered stress field is given over the whole parameter domain by the surrogate model 



(42) 
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where (5i)i^ii,h^} G H^"* ^re interpolation functionals, and are the only unknowns to be estimated in 
the "online" phase. 

Notice that gf''^{fi) satisfies the homogeneous equilibrium equations for any ii E V because the 

reduced basis stress fields (</'.)ie|i.s^] are linear combinations of fields that satisfy the homogeneous 

equilibrium equations. Then, by linearity, the recovered stress is indeed equilibrated in the finite 
element sense. 

Remark: In terms of implementation, the various stress fields introduced in this section are always 
stored at the quadrature points associated with the finite element space. Indeed, the integrals in which 
these fields are involved can be evaluated exactly using the finite element quadrature rule. This is due 
to the fact that these stress fields are all obtained by differentiation of finite element fields. 



For instance. To compute the POD operator in equation (40), we expand each of its component as 
follows: 



0~ ^ ^ ^ ~ 



git'''il4))--Dil^)--gi^''i&)dn 

+ I e(M^MD) : D{[4) Qipo) Dif^) : eSu'' {fJ^)) dQ 
Jn 



(43) 



This expansion is not used in practice. However, it shows that the components of the POD operator 
can be obtained exactly by applying the standard finite element quadrature rule. As a consequence, 
all the additive operations between stress fields, and the multiplication by scalars, can be performed 
at the quadrature point level. 

4.5 "Online" evaluation of the reduced order model for the stress 

In the "online" phase, a simple method to estimate the coefficients (5i(/i))ig|i_ft_^| of the recovered 
stress evaluated at a particular parameter values fi £ V would be to perform an a priori interpolation 
of these coefficients over the parameter domain. This could be done, for instance, by means of a moving 
least-squares technique or Kriging. 

A more advanced interpolation technique consists in computing a recovered stress d^{p) optimally 
in the sense of the maximisation of the efficiency of the error estimate i^"P(/i). In other words, for a 
particular /i S T', we look for a recovered stress field that is compatible with the surrogate model and 
is solution to the optimisation problem: 

gifj.)^ argmin 11^* -^h(/i)||c(/.) , (44) 

where the space of admissibility for the reduced stress field corresponding to parameter /i is defined 
by 5^(1};m) = {g* e S{n)\g* = f + ^''^^(m), V (a*).eii,f.,i e M"^- This problem of 

optimisation can be recast is the variational form 

Find |(/x) e 5' (17;/i) such that V^* G 5''° (17), 

K^i) : ciji) ■.g*dn+ [ I (u^iji)) ■.g*dn = o, ^"^^^ 
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where the constitutive relation has been used to obtain the right-hand side of the equation, and the 
space 5''0(f7) is defined by S'^^^in) = {g* e S{n)\g* = El^i^."*. V«)^eIl,,lJ e M"*}. 

Now, by recalling that y^{n) — +m'^''^(/^) and taking into account equation (34) and the fact 

that G Z^'^'"(fi), we obtain the following variational form for the determination of the recovered 

stress field: 

V^* e / : C'(m) ■.g'dn^ [ g {y}'•P{^J.)) : dQ . (46) 

Jn Jn 

Therefore, the expansion coefficients (Q?i(/i))ig|i are obtained by solving the linear system of 
equations 

K(a*)«(a*) = f'(^) (47) 
where the components of vector 5(/x) are the interpolation coefficients (ai{fi)) . , and 

V(*,j)e Il,n^F, k]U)= I 0. :C(M):0.df^ 

r r~' ~1 (48) 

Vj e II, nj, F;(Ai) = / e(M''P(M)) ■ t '^^ ■ 



"Offline" /"Online" split of the numerical complexity. Assembling linear system (47) can be 
done efficiently by a combination of "offline" and "online" computations, where all the operations 
with numerical complexity that depends on the dimension of the underlying finite element space are 
performed "offline". 

We will first assume that compliance C admits the natural affine form: 



VAie7',Vxel], C(x; m) = ^ ^(x) 7- (a^) • 

— r 

Then, the expression of operator K can be expanded as fohows: 



(49) 



k=\ 



where, V/c e Il,nel, V(z,i) e Il,n^f , K^^^^. 



(50) 



Similarly, by making use of the expansion of quantity u^'P over the parameter domain (equation ( 23 ) , 
we find that the right hand side of linear system ( 47 1 reads 



Va^gT', F (A.) = ^f;7fc(M) 



k = \ 



where, V fc e p , n^l , V j G p , , F^^^- = / ) : ^ dn 



(51) 



The "online" assembly operations, comprising the first line of (50) and (51 ) respectively, do not depend 
on the number of finite element unknowns. The remaining operations, namely the integrations specified 
by the second line of (50) and (51) respectively, are performed "offline". An expansion similar to (43) 
can be done to show that the offline terms can be integrated using the quadrature rule associated with 
the finite element space U^{il). 



15 



4.6 Computation of the upper bound 



Once coefficients {ai)i^fi^n^j corresponding to a particular parameter fi £ V have been computed, the 
"onhne" error estimate j^"p(/i) defined by 

= / (^^(m) - Km)) : Cip) ■■ ig'in) - gin)) ' (52) 
Jn 

can be evaluated. This can be done by a combination of "online" and "offline" computations that 
make use of the affine form of C, the affine expansions of a;'''" (equation (42 1 ) and a^'^ (equation ( [38| ) 

over the parameter domain, and the affine expansions of u'^''^ (equation (15)) and m'^'P (equation (23)). 
The technique is similar to the one deployed to assemble system (47) and will not be detailed for the 
sake of concision. 



4.7 Lower bound 

As will be shown in the results section of this paper, deriving an efficient lower bound for the "online" 
error associated to the reduced order modelling strategy can help control the efficiency of the upper 
bounding method. In particular, the choice of the dimension of the reduced space for the recovered 
stress remains to be made at this stage. 

A lower bound j/'°™(/i) for the error measure ||e' (/i)||Li(^) corresponding to any parameter ji £ V 

can be obtained by first constructing an enhanced surrogate for the displacement: 

VmGT', u^{p)~u^'{p)^u'''^{ii) + u^''\lj). (53) 

such that u2^'"(/i) belongs to a reduced space richer than W^'^iVL) [i.e. U''"{n) C U'^''"{fl) C 

U^''-\il)). In order to define this enriched space, we simply perform "offiine" a Snapshot-POD of higher 
order for the displacement field, which reads 

Z^2>-.0(r!):=span((^^),gIi,„..j) , (54) 
where Us > > and, consistently with the notations of section 



3.4 



(56) 



W^eln^ + l,nfl 0^ = £ ^h.0(^.) ^ (55) 

In the "online" phase, m^'^(m) is optimally obtained by a Galerkin projection of the governing 
equations in space U^'^'^{ft), which reads 

Find u^''°{ii) e U^'''"{n) such that Vu* e U^''°{n), 

The field u^'^ifJ,) is a hierarchically enhanced approximation of ^'^(/i). It is expected to be closer than 
u''(/i), in some sense, to the "truth" finite element solution u^ifi). The method to assemble and solve 
problem ( 56 1 using a combination of "offiine" and "online" operations is the same as the one detailed 
in section [3| the only difference being the dimensions of the respective reduced spaces. 

Theorem Defining an approximation of the error by ^(/i) '.— u^^in) — u''(/i) for any ^ £ V, we 
have the following lower bounding property 

\R(^i^i)■,^l)\ 

V^,eV, ^'°-(a.) := < lle'-(M)llD(,) , (57) 

\\e [t£)\\D(p.) - 

where the residual R is the parametrised linear form defined by 

y er,yu* £U{n), Riv;^l):^liu*]^i)-a{u\^l),u*;fJ.). (58) 
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Proof The proof is an extension of a similar technique used in the finite element context (see for 



instance [S])- The starting point is the weak form (13) of the parametrised problem of elasticity, from 
which we can obtain, for any fi Cz V, 

yu" eU^''"{n), aiy}'-°ifi)~u''°iii),u'';ii)^l{y^^ (59) 

Using the identities u'^ = u"^'" + vP and — y}^''^ + itP, we get the weak form of the equations governing 
the error e' : 

Vu^ e a{e'{^J),u*;lJ.)=R{u*:lI). (60) 

We can now substitute for u* the approximate error g '(^) e U^'"{n), which leads to the expression 

Vu* eU'^'^'in), a(e'-(M),r(Ai);A^)=i?(r(M);A^). (61) 

The exact error and its approximation e^(//) belong to the finite element space 14^''^{V,). The 

bilinear form a is an inner product for this particular space. We can therefore apply the Cauchy- 
Schwartz inequality to obtain 



^a{e^■{^i),e^{^i);^l)^a{^{^i),^■{^^) > |i?(r(/i);A^)l (62) 

and the announced result is immediate by making use of the definition of |j . O 

Again, an efficient combination of "offline" and "online" computations can be deployed such that 
the "online" numerical complexity associated with evaluating" z^'™ does not depend on the number of 
finite element unknowns. 

Notice that the lower bounding technique makes use of the Galerkin orthogonality. It is therefore 
restricted to the projection-based reduced order modelling scheme described in this paper, as opposed 
to the more general upper bounding methodology described previously. 

Summarising the upper and lower bounding results developed in the previous sections, we have, 
for any /i G P 

j.'°-(^)<||e^(M)||z,(^)<^^"P(A^). (63) 
4.8 Relative error measures 

Let us define the relative error bounds z/^p and by 

\u^'i[i)\\D{^) ' " IIw^'Xm)II^(m) 



V/iGP, ^"p-'"'(m):- „ 2r, f and ^i°w.rci(^) ^ ^ ^ (64) 



and the following properties are satisfied 

VfieV, ^'-■-'(/i)< „ ... ' <^"P'-^(a^). (65) 

A combination of "online" and "offline" computations can be deployed to compute the denominator 
of (64 1, by making use of the affine form of D and of the affine expansions of and w'^'P over the 



parameter domain. 
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5 Example: homogenisation of random composite materials 



5.1 Simple computational homogenisation framework for composite mate- 
rials 

We describe the example of two-dimensional homogenisation of a composite material for which we 
illustrate the proposed methodology. We consider an heterogeneous elastic material made of two phases 
with distinct elastic constants: circular inclusions and surrounding matrix (see figure[2|. The inclusions 
are distributed randomly in the material and have random diameters. The aim is to determine: 

• some quantities of interest of the local fields within the composite as a function of the far-field 
load applied to the material, parametrised by /i' € P', 

• the so-called effective elasticity tensor, which will be defined later on, as a function of some 
characteristics /x™ G of the material heterogeneities. 




Effective elasticity tensor Z)^ Effective strain 5 — {^2, ^J■3, 1^4} 




Figure 2: Schematic representation of the computational homogenisation framework for composite 
materials. The elastic constrast of the particulate composite is parametrised. 

Following classical approaches in homogenisation of random heterogeneous media (see for instance 
[17]), we consider that domain is a statistical volume element (SVE) of the composite material: 
a volume of simple shape (here square) of the composite material that is large enough to contain 
sufficient statistical information about the random distribution of inclusions. Over the SVE, the 
displacement field is additively split into a "coarse" contribution u and a fluctuation field u. Defining 
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fi = (^fi™^ fi^^^ , this reads: 

V/i e T', Va; e ri, u{x;fj,)—u{x]^) + u{x]^) (66) 
The "coarse" contribution u hes in a low-dimensional "smooth" space U{il.) of dimension 3 defined by 
U{n) = {u* e U{n) I u*{x) ^ e^{x-x),ye^ eR^ xR^ such that ^ i^^} , (67) 



where x is the barycenter of ft {i.e. Jq{x — x) dft ~ 0). Effective strain e'^ represents a far-field action 
applied to the material. The complementary fluctuation field w lies in space U{Q) of fields with small 
characteristic length of variation U{il) defined by: 

u{n) - [u* e uin) I {g{u*))n = o} , (68) 

where the averaging operator is defined { ■ ) fi ■= Jq ■ d^l. In order the fiuctuation field to be 
completely determined under the action of the coarse field, we need to define boundary conditions for 
the SVE problem. We choose to enforce Dirichlet boundary conditions on the boundary dfl = of 
the SVE, consistently with definition (68) and the so-called macrohomogeneity condition (see |27j): 

y lie'P,yx£ dn, w(x;/i) = u(z;^). (69) 

It is now possible, assuming that = to solve the parametrised boundary value problem {([7|), 



([8]), (69)} associated with the SVE to obtain the local displacement and stress fields as a function of 
the effective strain and the characteristics of the heterogeneities, which is called downscaling. 

Conversely, the upscaling operation consists in computing the effective elasticity tensor from 

the knowledge of the local fields, as a function of the distribution of the material heterogeneities. This 
effective tensor is defined by the relationship 



y^m ^ 7?m^v(^*,u*) e S^'^in) X (U{n)+U{n)'^ satisfying ([69| and ^, 



(70) 



The components of the effective elasticity tensor can be obtained by numerical testing (see for instance 
[40j ) . whereby one applies a range of elementary effective strains through the Dirichlet boundary 
conditions, solves the corresponding SVE boundary value problems and post-processes the resulting 
stress. 

Computing the local fields for arbitrary parameters of the heterogeneities and far-field action can 
require a tremendous numerical effort. We reduce this effort by deploying the Galerkin-POD described 
in section [3] for the parametrised SVE problem. 



5.2 Data and discretisation of the parametrised SVE problem 

The SVE boundary value problem {([t]), ([s]), ([69])} is parametrised by the material parameters ^™ and 
the load parameters /i'. Domain is the unit square in 2D, defined hy fl = [0, 1] x [0, 1], over which 
we discretise the elasticity problem with piecewise constant elasticity constants in a conforming way, 
using linear triangle elements (see figure [s]). 

We assume in this example that the material heterogeneities are only parametrised by the elastic 
contrast /i™ = iii — where iJ™^' is the Young's modulus of the matrix and i?'"'^ is the Young's 

modulus of the inclusions. The elastic contrast ranges from 0.1 (soft inclusions) to 10 (hard inclusions). 
The Poisson's ratios of both phases is set to v = 0.3. In this context, the affine representation of the 
Hooke's elasticity tensor over the parameter domain reads: 

yiier,yxen, d{x, ^i) ^ d'''^' + {^ii - 1) w'^^x) d""^' (71) 
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Figure 3: Finite element discretisation of the parametrised SVE problem. 



In the previous equations, function i?'"*^ is the indicator function of the inclusion phase. It is equal to 
1 for a point located in an inclusion and elsewhere. The elasticity tensor of the matrix phase is 

defined by equation ([2]), with — 1^ and v"^^^ — 0.3. The affine representation of the compliance 

tensor over the parameter domain becomes 

V^i e -P, Vx e n, C{x,fi) = C""^* + (^-^ - 1^ C""^* (72) 

where C™^* is the compliance tensor of the matrix phase. 

The independent components of the effective strain constitute parameters /i^ More precisely, we 
define ^\ = fj,2 = §^i, M2 = A*3 = £22 ^^'^ Ms = A'4 = £12- The affine representation of the Dirichlet 
boundary conditions are now 

V/^ e P, Vx e w(x,/i)==^^p 0) + 1 j A^s + q^M4^(^-s). (73) 

The parameter domain is restricted to the hypercube V ^ [0.1, 10] x [—1, 1] x [—1, 1] x [—1, 1]. 



5.3 Reduced order modelling and certification 

"Offline procedure" In order to deploy the Galerkin-POD methodology, we first sample the pa- 
rameter domain V. This is done by generating a 4-dimensional Quasi- Random sequence (the Sobol 
sequence implemented in MATLAB®) and retaining the first 20 points, which defines V (see figure 
|4]). The corresponding "truth" finite element displacement and stress fields are stored, for interpola- 
tion and error estimation purposes respectively. In our implementation, we store the nodal values of 
displacements, while we store the stresses at the quadrature points associated with U^{fl) {i.e.: at the 
centroids of the triangle elements). 
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Figure 4: Quasi Monte-Carlo Sampling of the parameter domain. Only three of the four dimensions 
are represented here. The blue points are the parameter values used to compute the snapshot. The red 
line is the ID manifold over which the "online" reduced order modelling interpolation and associated 
error estimation is validated. 



Next, we compute the set of lifting functions (V'j)ie|i.3l G (^^(^))'^- In order to do so, we first set 
fj,o — {I 0)'^ (we therefore make us of the bilinear form associated with an homogeneous material). 



and then solve the series of corresponding finite element problems (recall equation (22)) with Dirichlet 



boundary conditions corresponding to each of the term of affine expansion (73). We can then subtract 
the lifting u^'P(/i) from displacements w'^(/x) for all parameters fi belonging to the sampled parameter 
domain V. A Singular value decomposition is performed on the resulting bubble fields (they vanish 
over dfl = dft^). The first spatial modes of this proper orthogonal decomposition are then stored 
(we store their nodal values). 

The lifting for the stress needs not be performed as b{x) = 0, Va; € and dft^ = {}. A Snapshot 
POD in the sense of the energy norm associated to fio is performed on the set of sampled finite element 
stress, and the first spatial modes are stored (we store their values at the quadrature points). 

Last, we can pre-compute all the operators that allow for the "offline" / "online" split of the nu- 
merical complexity associated to interpolation and verification as described througout sections [3] and 

H 

"Online" procedure In the online phase, for a particular /x e P of interest, we need to evaluate the 
lifting functions w^'P(^) and (the latter is null in our case), and interpolate the complementary 

parts of the displacement and stress fields, and cr^'^ifi) respectively. The first set of these 



operations is done by direct evaluation (equation (22)), while the latter requires the solution of the 



projected problems for the displacement (equation (16)), for the hierarchically enriched displacement 



(equation (56)) and for the equilibrated stress in the finite element sense (equation (45)). 

Finally, the upper bound for the "online" error associated with the Galerkin-POD is obtained 
by comparison of the recovered stress field ct(/z) to the stress field obtained by applying the 

constitutive relation to the reduced approximation of the displacement (/i) , as explained in section 

tThe lower bound is obtained by comparing the hierarchically enriched displacement field M^'^(/i) to 
e reduced displacement u'Xa*)- 
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5.4 Numerical results 



We validate our error bounding methodology on a manifold of the parameter domain 7?°™' = g 
P I /Lt2 = 0, /i3 = and /i4 = 1}, which is illustrated in figure |4] and corresponds to applying a shear 
load to the SVE. Notice that none of the sampled parameter values actually belong to this manifold. 
The numerical results are given in figures [5] and [6] and are commented below. 

5.4.1 Effectivity of the error bounds 




Figure 5: Effectivity of the upper and lower bounds for the "online" error associated to the reduced 
order modelling approach as a function of the elastic contrast of the composite structure and of the 
expansion order of the two surrogates used to obtain the bounds. 

First, we show in figure [s] the influence of parameters and on the effectivity of the upper 
and lower bounds for the error introduced by the reduced order modelling approximation. 

For demonstration purposes, is arbitrarily set to 7. The left-hand side graph in figure [s] shows 

l|e'(Af)llD(A,) 

the evolution of the exact relative error in energy norm n , the upper bound iy"?''''^' for this 

error measure and the lower bound ^lo"''^^! as a function of the elastic contrast ^i. The upper bound 
is first obtained by using and expansion of order 5 for the stress surrogate {i.e. — 5). This number 
is then increased to 12 in an incremental manner. Similarly, the lower bound is first obtained using an 
expansion of order 8 for the enriched displacement surrogate {i.e. — 8), which is then increased to 
13 in an incremental manner. It clearly appears that both the lower and upper bounds consistently 
sharpen with increasing dimensions of the corresponding reduced spaces. 

The sharp decrease of the exact error around point (10 1)^ of the parameter domain can be 
explained by the fact that our Galerkin-POD model provides the exact finite element solution at this 
particular point. Indeed, the exact solution is obtained by direct evaluation of the lifting due 
to the choice that we made for /iq. The complementary part of the reduced solution vanishes 
at that point. 

The effectivity results stated previously are emphasised by the right-hand side graph in figure [Sj 
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The elastic contrast is fixed at a value of 1.6. The effectivity index of the two bounds, defined by 



\\e''{^J.)\\D^p) 

VA^eT', <; ,ow - (74) 



3 low 



are plotted as a function of the ratio between the the dimension of the corresponding reduced spaces, 
hif) and rv^ , and the dimension of the reduced space used to compute the initial approximation of the 
solution {i.e. = 7). Both bounds seem to converge quickly to an effectivity of 1. The lower bound 
is extremely effective even when using a very small number of additional POD modes. The upper 
bound needs slightly more computational effort to reach the expected numerical efficiency. However, 
the sharpness of the upper bound could surely be controlled in an adaptive manner using the available 
distance between the upper and the lower bounds. 



Remark: It should be noticed that in spite of the fact that the bounds seem to converge to the exact 
error with increasing dimensions of the corresponding reduced spaces, we have no reason to believe that 
this behaviour is asymptotically true. Indeed, the information that is available to train the surrogate 
models is obtained by sampling the parameter domain. Therefore, the exact finite element stress and 
exact finite element displacement cannot, in general, be obtained by simply increasing and . In 
fact, we should observe a stagnation of the effectivities when using reduced spaces of large dimensions 
to compute the bounds. To alleviate this potential issue, the only solution is to refine the sampling of 
the parameter domain in an adaptive manner. 



5.4.2 Convergence of the error bounds 

Next, we show how the proposed bounds converge with the error in the reduced order modelling 
approximation. To do so, we set = + 2 and = + 1. In this manner, the numerical 
parameters of the lower and upper bounds are constrained to follow the order of the surrogate model 

n ^ l|e''(/f)llD(A£) 

6 we show the convergence of the exact relative error n^r mi 



that we want to certify. In figure 



and the convergence of the upper and lower bounds z/"P'''°' and as a function of n^. As for the 

previous results, the plots are given for the manifold 7?'=™'. 

In the right-hand side bottom graph, we fix fii = 1.6 to ease the interpretation. These partial 
results show that the bounds have approximately the same convergence rate as the reduced order 
modelling approximation that we want to certify. 



6 Discussion and conclusion 

We have presented a novel and conceptually simple way to build lower and upper bounds for the error 
arising in POD-based reduced order models for parametrised elasticity problems. The upper bounding 
relies on the error in the constitutive relation, which requires to construct and evaluate a reduced order 
model for the stress field. This bound does not rely on the Galerkin orthogonality, which potentially 
allows us to apply it to other type of reduced order modelling techniques (Kriging, Moving Least- 
Squares approximations, response surface method, etc.). We have shown good numerical convergence 
and effectivity properties of both bounds, which are the basic requirements to perform adaptivity. 

In terms of implementation, we have shown that for each step of the bounding procedure, and 
adequate "offfine" / "online" split of the numerical complexity guarantees that the certification remains 
efficient numerically. However, the intrinsic nature of projection-based reduced order modelling makes 
the implementation complex, and highly intrusive. 
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Figure 6: Convergence of the upper and lower bounds of the reduced order modehing error as a fonction 
of the order of the POD expansion used to approximate the solution of the parametrised problem of 
computational homogenisation over the parameter domain. 



The bounding techniques have been applied to a basic example of multiscale modelling, for which 
the requirements in terms of cost reductions are still rather stringent. On the way, we have laid some 
foundations for future expansion of this work to more complex microstructures and micro/macro re- 
lationships. 

Our next step is to directly evaluate the surrogate error in terms of quantities of interest. In the 
case of elasticity, this is a rather straightforward step, as the goal-oriented error estimation techniques 
based on duality are very well-established. From there, we can imagine developing an adaptive reduced 
order modelling approach. To do that, efficient procedures need to be devised to control the dimension 
of the reduced space used to approximate the solution, the sharpness of the error bounds and the 
sampling of the parameter domain. 

In terms of homogenisation, the goal-oriented certification will require to devise two type of error 
bounds, one for the averaged information, controlled by a global quantity of interest, and one for the 
local information at the material level. The local quantities of interest are very often the maximum 
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value of some stress or strain components. Unfortunately, the corresponding dual problems are not 
reducible in the general case. An approximation will be necessary, at the probable cost of loosing the 
bounding property, or not bounding the right quantity. 

At last, we hope that the proposed approach complements existing work in the area of error 
estimation for reduced order modelling. A formal and numerical comparison of the methods available 
will be necessary. More importantly, we believe that this contribution will help progressing towards 
the certification of non-afiine and nonlinear parametrised boundary value problem, for instance by 
using the classical extensions of the the constitutive relation error to nonlinear problems. 
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